Step 3. Check analytical strategies (part 2)
# remove objects and memory
rm(list=ls());gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 835948 44.7 1627021 86.9 1178085 63
Vcells 1686107 12.9 8388608 64.0 3406660 26#remove images
while(!dev.cur())dev.off()
cat("\014")load(paste0(getwd(),"/_data/step2.RData"))#borrar caché
#system("fc-cache -f -v")
#check R version
if(Sys.info()["sysname"]=="Windows"){
if (getRversion() != "4.4.1") { stop("Requiere versión de R 4.4.1. Actual: ", getRversion()) }
}
if(Sys.info()["sysname"]=="Linux"){
if (getRversion() != "4.4.1") { stop("Requiere versión de R 4.4.1. Actual: ", getRversion()) }
}
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
# install.packages(c("dplyr", #for data
# "tidyr", #for data
# "lubridate", #for dates
# "openxlsx", #for excel files
# "rio", #for importing and exporting data
# "purrr", #for iterating in databases
# "devtools", #for external packages
# "DiagrammeR", #para visualizar DAG
# "dagitty",
# "ggdag",
# "ggplot2" #for graphics
# "kableExtra", #pretty tables
# "quarto", #for documents
# "geosphere" #for coordinates and classifying
# "geepack", #for regression
# "glmmTMB", #For GLMMs
# "DHARMa", #For residual diagnostics
# "car", #For hypothesis testing
# "brms", #Bayesian model
# "bayesplot",
# "loo",
# "Synth", #for synthetic control method
# "weathercan",#for weather data
# "sandwich", #cluster robust intervals
# "emmeans", #for predictions
# "gnm", #Conditional Poisson models
# "splines", #nonlinearity
# "geeM", #negative binomial and more flexible GEE models
# "PanelMatch", #Matching technique with panel data
# "scpi" #control sintético
# "nixtlar", #for time series analysis and prediction
# "CausalImpact" #for time series causal impact
# "forecast" #for time series analysis prediction and decomposition
# ))
library(dplyr); library(lubridate); library(tidyverse); library(openxlsx); library(rio); library(purrr); library(dagitty); library(ggdag); library(kableExtra); library(geosphere); library(geepack); library(lme4); library(glmmTMB); library(DHARMa); library(car); library(brms); library(bayesplot); library(loo); library(Synth); library(weathercan); library(sandwich); library(emmeans); library(gnm); library(splines); library(geeM); library(plm); library(PanelMatch); library(scpi); library(fect); library(nixtlar); library(CausalImpact); library(forecast)
Adjuntando el paquete: 'dplyr'The following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, union
Adjuntando el paquete: 'lubridate'The following objects are masked from 'package:base':
date, intersect, setdiff, unionWarning: package 'readr' was built under R version 4.4.3── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.2 ✔ tibble 3.2.1
✔ purrr 1.0.4 ✔ tidyr 1.3.1
✔ readr 2.1.5 ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Adjuntando el paquete: 'dagitty'
The following object is masked from 'package:rio':
convert
Adjuntando el paquete: 'ggdag'
The following object is masked from 'package:stats':
filter
Adjuntando el paquete: 'kableExtra'
The following object is masked from 'package:dplyr':
group_rows
Cargando paquete requerido: Matrix
Adjuntando el paquete: 'Matrix'
The following objects are masked from 'package:tidyr':
expand, pack, unpack
Adjuntando el paquete: 'lme4'
The following object is masked from 'package:rio':
factorize
This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
Cargando paquete requerido: carData
Adjuntando el paquete: 'car'
The following object is masked from 'package:purrr':
some
The following object is masked from 'package:dplyr':
recode
Cargando paquete requerido: Rcpp
Loading 'brms' package (version 2.22.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
Adjuntando el paquete: 'brms'
The following object is masked from 'package:glmmTMB':
lognormal
The following object is masked from 'package:lme4':
ngrps
The following object is masked from 'package:stats':
ar
This is bayesplot version 1.12.0
- Online documentation and vignettes at mc-stan.org/bayesplot
- bayesplot theme set to bayesplot::theme_default()
* Does _not_ affect other ggplot2 plots
* See ?bayesplot_theme_set for details on theme setting
Adjuntando el paquete: 'bayesplot'
The following object is masked from 'package:brms':
rhat
This is loo version 2.8.0
- Online documentation and vignettes at mc-stan.org/loo
- As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.
- Windows 10 users: loo may be very slow if 'mc.cores' is set in your .Rprofile file (see https://github.com/stan-dev/loo/issues/94).
##
## Synth Package: Implements Synthetic Control Methods.
## See https://web.stanford.edu/~jhain/synthpage.html for additional information.Warning: package 'weathercan' was built under R version 4.4.3As of v0.7.2, the `normals` column in `stations()` reflects whether or not there
are *any* normals available (not just the most recent).
Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'
Registered S3 method overwritten by 'lfe':
method from
nobs.felm broom
Adjuntando el paquete: 'plm'
The following objects are masked from 'package:dplyr':
between, lag, lead
Adjuntando el paquete: 'PanelMatch'
The following object is masked from 'package:tidyr':
extract
The following object is masked from 'package:stats':
weights
Registered S3 method overwritten by 'GGally':
method from
+.gg ggplot2
Cargando paquete requerido: bsts
Cargando paquete requerido: BoomSpikeSlab
Cargando paquete requerido: Boom
Adjuntando el paquete: 'Boom'
The following objects are masked from 'package:brms':
ddirichlet, rdirichlet
The following object is masked from 'package:stats':
rWishart
Adjuntando el paquete: 'BoomSpikeSlab'
The following object is masked from 'package:stats':
knots
Cargando paquete requerido: zoo
Adjuntando el paquete: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
Cargando paquete requerido: xts
######################### Warning from 'xts' package ##########################
# #
# The dplyr lag() function breaks how base R's lag() function is supposed to #
# work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
# source() into this session won't work correctly. #
# #
# Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
# conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
# dplyr from breaking base R's lag() function. #
# #
# Code in packages is not affected. It's protected by R's namespace mechanism #
# Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
# #
###############################################################################
Adjuntando el paquete: 'xts'
The following objects are masked from 'package:dplyr':
first, last
Adjuntando el paquete: 'bsts'
The following object is masked from 'package:BoomSpikeSlab':
SuggestBurn
Registered S3 method overwritten by 'quantmod':
method from
as.zoo.data.frame zoo
Adjuntando el paquete: 'forecast'
The following object is masked from 'package:brms':
ma#special repository indicated or the package
if(!require(weathercan)){
install.packages("weathercan",
repos = c("https://ropensci.r-universe.dev", "https://cloud.r-project.org")); library(weathercan)
}
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
if(!require(bpmn)){devtools::install_github("bergant/bpmn")}Cargando paquete requerido: bpmn#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
sum_dates <- function(x){
cbind.data.frame(
min= as.Date(min(unclass(as.Date(x)), na.rm=T), origin = "1970-01-01"),
p001= as.Date(quantile(unclass(as.Date(x)), .001, na.rm=T), origin = "1970-01-01"),
p005= as.Date(quantile(unclass(as.Date(x)), .005, na.rm=T), origin = "1970-01-01"),
p025= as.Date(quantile(unclass(as.Date(x)), .025, na.rm=T), origin = "1970-01-01"),
p25= as.Date(quantile(unclass(as.Date(x)), .25, na.rm=T), origin = "1970-01-01"),
p50= as.Date(quantile(unclass(as.Date(x)), .5, na.rm=T), origin = "1970-01-01"),
p75= as.Date(quantile(unclass(as.Date(x)), .75, na.rm=T), origin = "1970-01-01"),
p975= as.Date(quantile(unclass(as.Date(x)), .975, na.rm=T), origin = "1970-01-01"),
p995= as.Date(quantile(unclass(as.Date(x)), .995, na.rm=T), origin = "1970-01-01"),
p999= as.Date(quantile(unclass(as.Date(x)), .999, na.rm=T), origin = "1970-01-01"),
max= as.Date(max(unclass(as.Date(x)), na.rm=T), origin = "1970-01-01")
)
}
smd_bin <- function(x,y){
z <- x*(1-x)
t <- y*(1-y)
k <- sum(z,t)
l <- k/2
return((x-y)/sqrt(l))
}
theme_custom_sjplot2 <- function(base_size = 12, base_family = "") {
theme_minimal(base_size = base_size, base_family = base_family) +
theme(
# Text elements
text = element_text(size = base_size, family = base_family),
plot.title = element_text(face = "bold", hjust = 0.5, size = base_size * 1.2),
plot.subtitle = element_text(hjust = 0.5, margin = margin(b = 10)),
axis.title = element_text(size = base_size, face = "bold"),
axis.text = element_text(size = base_size * 0.8),
axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.5),
axis.text.y = element_text(angle = 0, hjust = 1, vjust = 0.5),
axis.title.x = element_text(margin = margin(t = 10)),
axis.title.y = element_text(margin = margin(r = 10)),
# Plot layout
plot.margin = margin(t = 20, r = 20, b = 20, l = 20),
panel.grid.major = element_line(color = "grey80"),
panel.grid.minor = element_blank(),
legend.position = "right",
legend.text = element_text(size = base_size * 0.8),
legend.title = element_text(size = base_size, face = "bold"),
legend.background = element_rect(fill = "white", colour = NA),
legend.box.background = element_rect(colour = "grey80", linetype = "solid"),
legend.key = element_rect(fill = "white", colour = "white")
)
}
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
num_cores <- parallel::detectCores() -1
data.table::setDTthreads(threads = num_cores)#restore_after_fork = NULL, throttle = NULL)
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#CONFIG #######################################################################
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
options(scipen=2) #display numbers rather scientific number
nixtlar::nixtla_client_setup(api_key = readLines(paste0(gsub("f1/","", getwd() ),"/key.txt"))[[2]])Warning in readLines(paste0(gsub("f1/", "", getwd()), "/key.txt")): incomplete
final line found on 'H:/Mi unidad/PERSONAL
ANDRES/UCH_salud_publica/pasantia/f1/key.txt'API key has been set for the current session.#https://cran.r-project.org/web/packages/nixtlar/vignettes/prediction-intervals.html
#https://cran.r-project.org/web/packages/nixtlar/vignettes/exogenous-variables.html
#https://cran.r-project.org/web/packages/nixtlar/vignettes/cross-validation.html
#https://cran.r-project.org/web/packages/nixtlar/vignettes/anomaly-detection.html
collisions_weather_corr_rect_synth_tr$Y_lag10 <-
collisions_weather_corr_rect_synth_tr|>
mutate(Y_lag10 = dplyr::lag(rate_veh, 10, order_by = yday_corr))|>
mutate(Y_lag10 = dplyr::lag(rate_veh, 10, order_by = yday_corr))|>
ungroup()|>
pull(Y_lag10)
collisions_weather_corr_rect_synth_tr$Y_lag10 <-
collisions_weather_corr_rect_synth_tr|>
mutate(
third_value = nth(na.omit(Y_lag10), 1),
Y_lag10 = ifelse(is.na(Y_lag10), third_value, Y_lag10)
)|>
dplyr::select(Y_lag10)|>
ungroup()|>
pull(Y_lag10)
df35 <- collisions_weather_corr_rect_synth_tr %>%
arrange(tr_contr_sens, year.x, yday_corr) %>% # 1,2,3… days relative to race
filter(yday_corr<36) |>
transmute(
unique_id = "treated", # exactly this name
ds = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
y = rate_veh,
mean_median_lag_2_prec_median_imp = mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin= mean_median_total_precip_median_lin,
mean_min_temp_mean_lin= mean_min_temp_mean_lin,
mean_max_temp_mean_lin= mean_max_temp_mean_lin,
Y_lag10= Y_lag10
)
df35_post <- collisions_weather_corr_rect_synth_tr %>%
arrange(tr_contr_sens, year.x, yday_corr) %>% # 1,2,3… days relative to race
filter(yday_corr>=36, yday_corr<=max_time) |> # take the first 30
transmute(
unique_id = "treated", # exactly this name
ds = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
y = rate_veh,
mean_median_lag_2_prec_median_imp = mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin= mean_median_total_precip_median_lin,
mean_min_temp_mean_lin= mean_min_temp_mean_lin,
mean_max_temp_mean_lin= mean_max_temp_mean_lin,
Y_lag10= Y_lag10
)
df33 <- collisions_weather_corr_rect_synth_tr %>%
arrange(tr_contr_sens, year.x, yday_corr) %>% # 1,2,3… days relative to race
filter(yday_corr<34) |>
transmute(
unique_id = "treated", # exactly this name
ds = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
y = rate_veh,
mean_median_lag_2_prec_median_imp = mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin= mean_median_total_precip_median_lin,
mean_min_temp_mean_lin= mean_min_temp_mean_lin,
mean_max_temp_mean_lin= mean_max_temp_mean_lin,
Y_lag10= Y_lag10
)
df33_post <- collisions_weather_corr_rect_synth_tr %>%
arrange(tr_contr_sens, year.x, yday_corr) %>% # 1,2,3… days relative to race
filter(yday_corr>=34, yday_corr<=max_time) |> # take the first 30
transmute(
unique_id = "treated", # exactly this name
ds = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
y = rate_veh,
mean_median_lag_2_prec_median_imp = mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin= mean_median_total_precip_median_lin,
mean_min_temp_mean_lin= mean_min_temp_mean_lin,
mean_max_temp_mean_lin= mean_max_temp_mean_lin,
Y_lag10= Y_lag10
)
fc <- nixtla_client_forecast(
df35,
h = 7,
level = c(95),
freq = "D", # <— this fixes the infer_frequency error
finetune_steps = 10, #Number of steps used to finetune 'TimeGPT' in the new data.
finetune_loss ="rmse"
)Using historical exogenous features: mean_median_lag_2_prec_median_imp, mean_median_total_precip_median_lin, mean_min_temp_mean_lin, mean_max_temp_mean_lin, Y_lag10fc33 <- nixtla_client_forecast(
rbind(mutate(df33[1:2,], ds=ds-3), df33),
h = 7,
level = c(95),
freq = "D", # <— this fixes the infer_frequency error
finetune_steps = 10, #Number of steps used to finetune 'TimeGPT' in the new data.
finetune_loss ="rmse"
)Using historical exogenous features: mean_median_lag_2_prec_median_imp, mean_median_total_precip_median_lin, mean_min_temp_mean_lin, mean_max_temp_mean_lin, Y_lag10p <- nixtla_client_plot(rbind(df35, df35_post), fc,max_insample_length = 100)
p33 <- nixtla_client_plot(rbind(mutate(df33[1:2,], ds=ds-3), df33, df33_post), fc33,
max_insample_length = 100)
# 1) Banda de CI (layer 1)
p$layers[[1]]$geom_params$fill <- "#FFC10740" # amarillo semitransparente
p$layers[[1]]$geom_params$alpha <- 0.4
# 2) Línea histórica y pronóstico (layer 2)
# p$layers[[2]]$aes_params$colour <- "steelblue"
p$layers[[2]]$aes_params$linewidth <- 1.0
p33$layers[[1]]$aes_params$fill <- "#eb8f8f" # verde semitransparente
p33$layers[[1]]$aes_params$alpha <- 0.5
# 1) Banda de CI (layer 1)
p33$layers[[1]]$geom_params$fill <- "#FFC10740" # amarillo semitransparente
p33$layers[[1]]$geom_params$alpha <- 0.4
# 2) Línea histórica y pronóstico (layer 2)
# p$layers[[2]]$aes_params$colour <- "steelblue"
p33$layers[[2]]$aes_params$linewidth <- 1.0
p33$layers[[1]]$aes_params$fill <- "#eb8f8f" # verde semitransparente
p33$layers[[1]]$aes_params$alpha <- 0.5
t0 <- min(p$data$ds)
t0_33 <- min(p33$data$ds)
t0_posix <- min(p$data$ds)
t0_33_posix <- min(p33$data$ds)
xint35 <- t0_posix + days(34)
t0 <- as.Date(min(p$data$ds))
t0_33 <- as.Date(min(p33$data$ds))
# 3) Aplica scale_x_datetime() para transformar POSIXct → numérico
p +
scale_color_manual(
values = c(
# y = "#FFD600", # amarillo intenso
TimeGPT = "red" # violeta (DarkOrchid / BlueViolet)
),
labels= c( TimeGPT = "Prediction (95% IC in gray area)" )
) +
scale_x_datetime(
name = "Days of follow-up (35 days before race,\n7 days after the race)",
limits = c(min(p$data$ds), max(p$data$ds)),
breaks = seq(min(p$data$ds), max(p$data$ds), by = "7 days"),
labels = function(x) {
# convierto POSIXct a Date y resto t0
as.integer(as.Date(x) - t0 + 1)
}
) +
theme_classic(base_size = 13) +
theme(
strip.text = element_blank(), # quita el label de faceta
strip.background = element_blank(),
legend.position= "bottom"
)+
geom_vline(
xintercept = xint35,
linetype = "dashed",
color = "black",
linewidth = 1
) +
labs(y="Highspeed collisions per 1M vehicles")Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
p33+
scale_color_manual(
values = c(
# y = "#FFD600", # amarillo intenso
TimeGPT = "red" # violeta (DarkOrchid / BlueViolet)
),
labels= c( TimeGPT = "Prediction (95% IC in gray area)" )
) +
scale_x_datetime(
name = "Days of follow-up (33 days before race,\n7 days after the race)",
limits = c(min(p33$data$ds), max(p33$data$ds)),
breaks = seq(min(p33$data$ds), max(p33$data$ds), by = "7 days"),
labels = function(x) {
# convierto POSIXct a Date y resto t0
as.integer(as.Date(x) - t0_33 + 1)
}
) +
theme_classic(base_size = 13) +
theme(
strip.text = element_blank(), # quita el label de faceta
strip.background = element_blank(),
legend.position= "bottom"
)+
geom_vline(
xintercept = xint35-days(2),
linetype = "dashed",
color = "black",
linewidth = 1
) +
labs(y="Highspeed collisions per 1M vehicles")Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
## ── Lags ──
collisions_weather_corr_rect_quebec_synth_tr$Y_lag10 <-
collisions_weather_corr_rect_quebec_synth_tr |>
arrange(tr_contr, year.x, yday_corr) |>
mutate(
Y_lag10 = dplyr::lag(rate_veh, 10, order_by = yday_corr),
Y_lag10 = ifelse(is.na(Y_lag10),
first(na.omit(Y_lag10)),
Y_lag10)
) |>
ungroup() |>
pull(Y_lag10)
## ── Data frames (35 y 35 + 7) ──
df35_qc <- collisions_weather_corr_rect_quebec_synth_tr %>%
arrange(tr_contr, year.x, yday_corr) %>%
filter(yday_corr < 36) %>%
transmute(
unique_id = "treated",
ds = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
y = rate_veh,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin,
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
Y_lag10
)
df35_post_qc <- collisions_weather_corr_rect_quebec_synth_tr %>%
arrange(tr_contr, year.x, yday_corr) %>%
filter(yday_corr >= 36, yday_corr <= max_time) %>%
transmute(
unique_id = "treated",
ds = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
y = rate_veh,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin,
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
Y_lag10
)
## ── TimeGPT ──
fc_qc <- nixtla_client_forecast(
df35_qc,
h = 7,
level = c(95),
freq = "D",
finetune_steps = 10,
finetune_loss = "rmse"
)Using historical exogenous features: mean_median_lag_2_prec_median_imp, mean_median_total_precip_median_lin, mean_min_temp_mean_lin, mean_max_temp_mean_lin, Y_lag10p_qc <- nixtla_client_plot(
rbind(df35_qc, df35_post_qc),
fc_qc,
max_insample_length = 100
)
## Personalización idéntica a tu original
p_qc$layers[[1]]$geom_params$fill <- "#FFC10740"
p_qc$layers[[1]]$geom_params$alpha <- 0.4
p_qc$layers[[2]]$aes_params$linewidth <- 1.0
t0_qc <- as.Date(min(p_qc$data$ds))
xint35_qc <- min(p_qc$data$ds) + lubridate::days(34)
p_qc +
scale_color_manual(
values = c(TimeGPT = "red"),
labels = c(TimeGPT = "Prediction (95% IC in gray area)")
) +
scale_x_datetime(
name = "Days of follow-up (35 before, 7 after)",
limits = c(min(p_qc$data$ds), max(p_qc$data$ds)),
breaks = seq(min(p_qc$data$ds), max(p_qc$data$ds), by = "7 days"),
labels = function(x) as.integer(as.Date(x) - t0_qc + 1)
) +
theme_classic(base_size = 13) +
theme(
strip.text = element_blank(),
strip.background = element_blank(),
legend.position = "bottom"
) +
geom_vline(
xintercept = xint35_qc,
linetype = "dashed",
color = "black",
linewidth = 1
) +
labs(y = "High-speed collisions per 1 M vehicles")Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
## ── Y_lag10 en un solo pipe ───────────────────────────────────────────────────
collisions_weather_corr_rect_quebec_synth_tr$Y_lag10 <-
collisions_weather_corr_rect_quebec_synth_tr |>
arrange(tr_contr, year.x, yday_corr) |>
mutate(
Y_lag10 = dplyr::lag(rate_veh, 10, order_by = yday_corr),
Y_lag10 = ifelse(is.na(Y_lag10),
first(na.omit(Y_lag10)),
Y_lag10)
) |>
ungroup() |>
pull(Y_lag10)
## ── Conjuntos 33 días pre-intervención y 7 post ───────────────────────────────
df33_qc <- collisions_weather_corr_rect_quebec_synth_tr %>%
arrange(tr_contr, year.x, yday_corr) %>%
filter(yday_corr < 34) %>%
transmute(
unique_id = "treated",
ds = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
y = rate_veh,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin,
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
Y_lag10
)
df33_post_qc <- collisions_weather_corr_rect_quebec_synth_tr %>%
arrange(tr_contr, year.x, yday_corr) %>%
filter(yday_corr >= 34, yday_corr <= max_time) %>%
transmute(
unique_id = "treated",
ds = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
y = rate_veh,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin,
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
Y_lag10
)
## ── TimeGPT ───────────────────────────────────────────────────────────────────
fc33_qc <- nixtla_client_forecast(
rbind(mutate(df33_qc[1:2,], ds = ds - 3), df33_qc),
h = 7,
level = c(95),
freq = "D",
finetune_steps = 10,
finetune_loss = "rmse"
)Using historical exogenous features: mean_median_lag_2_prec_median_imp, mean_median_total_precip_median_lin, mean_min_temp_mean_lin, mean_max_temp_mean_lin, Y_lag10p33_qc <- nixtla_client_plot(
rbind(mutate(df33_qc[1:2,], ds = ds - 3),
df33_qc,
df33_post_qc),
fc33_qc,
max_insample_length = 100
)
## ── Estética idéntica a tu original ───────────────────────────────────────────
p33_qc$layers[[1]]$geom_params$fill <- "#FFC10740"
p33_qc$layers[[1]]$geom_params$alpha <- 0.4
p33_qc$layers[[2]]$aes_params$linewidth <- 1.0
t0_33_qc <- as.Date(min(p33_qc$data$ds))
xint33_qc <- min(p33_qc$data$ds) + lubridate::days(32)
p33_qc +
scale_color_manual(
values = c(TimeGPT = "red"),
labels = c(TimeGPT = "Prediction (95% IC in gray area)")
) +
scale_x_datetime(
name = "Days of follow-up (33 before, 7 after)",
limits = c(min(p33_qc$data$ds), max(p33_qc$data$ds)),
breaks = seq(min(p33_qc$data$ds), max(p33_qc$data$ds), by = "7 days"),
labels = function(x) as.integer(as.Date(x) - t0_33_qc + 1)
) +
theme_classic(base_size = 13) +
theme(
strip.text = element_blank(),
strip.background = element_blank(),
legend.position = "bottom"
) +
geom_vline(
xintercept = xint33_qc,
linetype = "dashed",
color = "black",
linewidth = 1
) +
labs(y = "High-speed collisions per 1 M vehicles")Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
We used the CausalImpact package, useful for time series data, inferring a counterfactual post race period based on pre-intervention components1.
# Model 2
ss2d <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2d <- AddLocalLevel(ss2d,c(as.numeric(df35$y), rep(NA,7))
) #
# Add weekly seasonal
ss2d <- AddSeasonal(ss2d,c(as.numeric(df35$y), rep(NA,7)),
nseasons=7, season.duration = 1) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2d <- AddAutoAr(ss2d, y = c(as.numeric(df35$y), rep(NA,7)), lags = 10) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
# Add external regressors explicitly (e.g., temperatures, precipitation)
covariates <- cbind(
# y = as.numeric(df35$y),
mean_min_temp_mean_lin = df35$mean_min_temp_mean_lin,
mean_max_temp_mean_lin = df35$mean_max_temp_mean_lin,
mean_median_lag_2_prec_median_imp = df35$mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin = df35$mean_median_total_precip_median_lin,
Y_lag10 = df35$Y_lag10
)
# Ensure names explicitly
colnames(covariates) <- c(#"y",
"mean_min_temp_mean_lin",
"mean_max_temp_mean_lin",
"mean_median_lag_2_prec_median_imp",
"mean_median_total_precip_median_lin",
"Y_lag10")
covariates<- as.data.frame(covariates)
covariates_future <- as.matrix(df35_post[, c("mean_min_temp_mean_lin",
"mean_max_temp_mean_lin",
"mean_median_lag_2_prec_median_imp",
"mean_median_total_precip_median_lin",
"Y_lag10")])
covariates_full <- rbind(covariates, covariates_future)
# Include regression terms
ss2d <- AddDynamicRegression(ss2d, covariates_full)
model2d_ratio <- bsts(c(as.numeric(df35$y), rep(NA,7)),
state.specification = ss2d, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
family ="student", #A Bayesian Analysis of Time-Series Event Count Data
niter = 3e4,
#burn = 2000, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
seed= 2125)
## Since we have predictor variables in the model, we need to explicitly
# make their coefficients time-varying using AddDynamicRegression().
impact3d_ratio_resp <- CausalImpact(bsts.model = model2d_ratio,model.args = list(prior.level.sd=.1, dynamic.regression=T),
post.period.response = df35_post$y)
plot(impact3d_ratio_resp, "original") 
summary(impact3d_ratio_resp)
plot(impact3d_ratio_resp)
=-=-=-=-= Iteration 0 Tue Jul 8 00:27:32 2025
=-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul 8 00:27:35 2025
=-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul 8 00:27:38 2025
=-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul 8 00:27:41 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul 8 00:27:44 2025
=-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul 8 00:27:47 2025
=-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul 8 00:27:50 2025
=-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul 8 00:27:53 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul 8 00:27:56 2025
=-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul 8 00:27:59 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 1.5 10.5
Prediction (s.d.) 1.6 (0.49) 11.5 (3.41)
95% CI [0.7, 2.6] [4.9, 18.3]
Absolute effect (s.d.) -0.14 (0.49) -0.97 (3.41)
95% CI [-1.1, 0.81] [-7.8, 5.65]
Relative effect (s.d.) -1.3% (766%) -1.3% (766%)
95% CI [-43%, 113%] [-43%, 113%]
Posterior tail-area probability p: 0.38543
Posterior prob. of a causal effect: 61%
For more details, type: summary(impact, "report")We included other terms (local linear trend).
library(bsts)
y35 <- c(as.numeric(df35$y), rep(NA, 7)) # respuesta escalada por 100 millones
ss2e <- list()
ss2e <- AddLocalLinearTrend(ss2e, y35)
ss2e <- AddSeasonal(ss2e, y35, nseasons = 7, season.duration = 1)
ss2e <- AddAutoAr(ss2e, y35, lags = 1)
ss2e <- AddDynamicRegression(ss2e, covariates_full)
set.seed(2125)
model35 <- bsts(y35, state.specification = ss2e, seed=2125,
niter = 3e4, family = "student")
#plot(model, "components", burn = 2000)
impact3d_ratio_resp2 <- CausalImpact(bsts.model = model35,model.args = list(prior.level.sd=.1, dynamic.regression=T),
post.period.response = df35_post$y)
#plot(impact3d, "original")
summary(impact3d_ratio_resp2)
plot(impact3d_ratio_resp2)
=-=-=-=-= Iteration 0 Tue Jul 8 00:28:10 2025
=-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul 8 00:28:13 2025
=-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul 8 00:28:15 2025
=-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul 8 00:28:18 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul 8 00:28:21 2025
=-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul 8 00:28:23 2025
=-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul 8 00:28:26 2025
=-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul 8 00:28:28 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul 8 00:28:30 2025
=-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul 8 00:28:33 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 1.5 10.5
Prediction (s.d.) 1.7 (0.53) 11.8 (3.70)
95% CI [0.67, 2.8] [4.69, 19.5]
Absolute effect (s.d.) -0.19 (0.53) -1.30 (3.70)
95% CI [-1.3, 0.83] [-8.9, 5.83]
Relative effect (s.d.) 12% (2224%) 12% (2224%)
95% CI [-46%, 119%] [-46%, 119%]
Posterior tail-area probability p: 0.35602
Posterior prob. of a causal effect: 64%
For more details, type: summary(impact, "report")# Model 2
ss2a <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2a <- AddLocalLevel(ss2a,c(as.numeric(df33$y), rep(NA,7))
) #
# Add weekly seasonal
ss2a <- AddSeasonal(ss2a,c(as.numeric(df33$y), rep(NA,7)),
nseasons=7, season.duration = 1) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2a <- AddAutoAr(ss2a, y = c(as.numeric(df33$y), rep(NA,7)), lags = 10) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
# Add external regressors explicitly (e.g., temperatures, precipitation)
covariates33 <- cbind(
# y = as.numeric(df35$y),
mean_min_temp_mean_lin = df33$mean_min_temp_mean_lin,
mean_max_temp_mean_lin = df33$mean_max_temp_mean_lin,
mean_median_lag_2_prec_median_imp = df33$mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin= df33$mean_median_total_precip_median_lin,
Y_lag10 = df33$Y_lag10
)
# Ensure names explicitly
colnames(covariates33) <- c(#"y",
"mean_min_temp_mean_lin",
"mean_max_temp_mean_lin",
"mean_median_lag_2_prec_median_imp",
"mean_median_total_precip_median_lin",
"Y_lag10")
covariates33<- as.data.frame(covariates33)
covariates_future33 <- as.matrix(df33_post[1:7, c("mean_min_temp_mean_lin",
"mean_max_temp_mean_lin",
"mean_median_lag_2_prec_median_imp",
"mean_median_total_precip_median_lin",
"Y_lag10")])
covariates_full33 <- rbind(covariates33, covariates_future33)
# Include regression terms
ss2a <- AddDynamicRegression(ss2a, covariates_full33)
model2a_ratio <- bsts(c(as.numeric(df33$y), rep(NA,7)),
state.specification = ss2a, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
family ="student", #A Bayesian Analysis of Time-Series Event Count Data
niter = 3e4,
#burn = 2000, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
seed= 2125)
## Since we have predictor variables in the model, we need to explicitly
# make their coefficients time-varying using AddDynamicRegression().
impact3a_ratio_resp <- CausalImpact(bsts.model = model2a_ratio,model.args = list(prior.level.sd=.1, dynamic.regression=T),
post.period.response = df33_post$y[1:7])
plot(impact3a_ratio_resp, "original") 
summary(impact3a_ratio_resp)
plot(impact3a_ratio_resp)
=-=-=-=-= Iteration 0 Tue Jul 8 00:28:40 2025
=-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul 8 00:28:43 2025
=-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul 8 00:28:46 2025
=-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul 8 00:28:49 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul 8 00:28:52 2025
=-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul 8 00:28:55 2025
=-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul 8 00:28:58 2025
=-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul 8 00:29:01 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul 8 00:29:04 2025
=-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul 8 00:29:07 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 1.5 10.4
Prediction (s.d.) 1.6 (0.43) 11.0 (2.98)
95% CI [0.74, 2.4] [5.19, 16.9]
Absolute effect (s.d.) -0.086 (0.43) -0.599 (2.98)
95% CI [-0.94, 0.74] [-6.55, 5.17]
Relative effect (s.d.) 2.2% (258%) 2.2% (258%)
95% CI [-39%, 99%] [-39%, 99%]
Posterior tail-area probability p: 0.4221
Posterior prob. of a causal effect: 58%
For more details, type: summary(impact, "report")library(bsts)
y33b <- c(as.numeric(df33$y), rep(NA, 7)) # respuesta escalada por 100 millones
ss2b <- list()
ss2b <- AddLocalLinearTrend(ss2b, y33b)
ss2b <- AddSeasonal(ss2b, y33b, nseasons = 7, season.duration = 1)
ss2b <- AddAutoAr(ss2b, y33b, lags = 1)
ss2b <- AddDynamicRegression(ss2b, covariates_full33)
set.seed(2125)
model33b <- bsts(y33b, state.specification = ss2b, seed=2125,
niter = 3e4, family = "student")
#plot(model, "components", burn = 2000)
impact3d_ratio_resp2 <- CausalImpact(bsts.model = model33b,model.args = list(prior.level.sd=.1, dynamic.regression=T),
post.period.response = df33_post$y[1:7])
#plot(impact3d, "original")
summary(impact3d_ratio_resp2)
plot(impact3d_ratio_resp2)
=-=-=-=-= Iteration 0 Tue Jul 8 00:29:17 2025
=-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul 8 00:29:20 2025
=-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul 8 00:29:22 2025
=-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul 8 00:29:25 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul 8 00:29:27 2025
=-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul 8 00:29:30 2025
=-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul 8 00:29:32 2025
=-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul 8 00:29:34 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul 8 00:29:37 2025
=-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul 8 00:29:40 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 1.5 10.4
Prediction (s.d.) 1.6 (0.47) 11.4 (3.27)
95% CI [0.74, 2.6] [5.17, 18.1]
Absolute effect (s.d.) -0.15 (0.47) -1.06 (3.27)
95% CI [-1.1, 0.74] [-7.7, 5.20]
Relative effect (s.d.) -6.3% (674%) -6.3% (674%)
95% CI [-43%, 99%] [-43%, 99%]
Posterior tail-area probability p: 0.37243
Posterior prob. of a causal effect: 63%
For more details, type: summary(impact, "report")# Model 2
ss2ab <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2a <- AddLocalLevel(ss2ab,c(as.numeric(df33$y), rep(NA,3))
) #
# Add weekly seasonal
ss2ab <- AddSeasonal(ss2ab,c(as.numeric(df33$y), rep(NA,3)),
nseasons=7, season.duration = 1) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2ab <- AddAutoAr(ss2ab, y = c(as.numeric(df33$y), rep(NA,3)), lags = 10) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
covariates_future33b <- as.matrix(df33_post[1:3, c("mean_min_temp_mean_lin",
"mean_max_temp_mean_lin",
"mean_median_lag_2_prec_median_imp",
"mean_median_total_precip_median_lin",
"Y_lag10")])
covariates_full33b <- rbind(covariates33, covariates_future33b)
# Include regression terms
ss2ab <- AddDynamicRegression(ss2ab, covariates_full33b)
model2ab_ratio <- bsts(c(as.numeric(df33$y), rep(NA,3)),
state.specification = ss2ab, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
family ="student", #A Bayesian Analysis of Time-Series Event Count Data
niter = 3e4,
#burn = 2000, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
seed= 2125)
## Since we have predictor variables in the model, we need to explicitly
# make their coefficients time-varying using AddDynamicRegression().
impact3ab_ratio_resp <- CausalImpact(bsts.model = model2ab_ratio,model.args = list(prior.level.sd=.1, dynamic.regression=T),
post.period.response = df33_post$y[1:3])
plot(impact3ab_ratio_resp, "original") 
summary(impact3ab_ratio_resp)
plot(impact3ab_ratio_resp)
=-=-=-=-= Iteration 0 Tue Jul 8 00:29:47 2025
=-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul 8 00:29:49 2025
=-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul 8 00:29:52 2025
=-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul 8 00:29:54 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul 8 00:29:57 2025
=-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul 8 00:29:59 2025
=-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul 8 00:30:01 2025
=-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul 8 00:30:04 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul 8 00:30:06 2025
=-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul 8 00:30:09 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 1.6 4.8
Prediction (s.d.) 1.7 (0.35) 5.1 (1.04)
95% CI [1.0, 2.4] [3.1, 7.2]
Absolute effect (s.d.) -0.11 (0.35) -0.34 (1.04)
95% CI [-0.81, 0.56] [-2.42, 1.68]
Relative effect (s.d.) -2.1% (26%) -2.1% (26%)
95% CI [-34%, 54%] [-34%, 54%]
Posterior tail-area probability p: 0.37121
Posterior prob. of a causal effect: 63%
For more details, type: summary(impact, "report")library(bsts)
y33b <- c(as.numeric(df33$y), rep(NA, 3)) # respuesta escalada por 100 millones
ss2b <- list()
ss2b <- AddLocalLinearTrend(ss2b, y33b)
ss2b <- AddSeasonal(ss2b, y33b, nseasons = 7, season.duration = 1)
ss2b <- AddAutoAr(ss2b, y33b, lags = 1)
ss2b <- AddDynamicRegression(ss2b, covariates_full33b)
set.seed(2125)
model33b <- bsts(y33b, state.specification = ss2b, seed=2125,
niter = 3e4, family = "student")
#plot(model, "components", burn = 2000)
impact3db_ratio_resp2 <- CausalImpact(bsts.model = model33b, model.args = list(prior.level.sd=.1, dynamic.regression=T),
post.period.response = df33_post$y[1:3])
#plot(impact3d, "original")
summary(impact3db_ratio_resp2)
plot(impact3db_ratio_resp2)
=-=-=-=-= Iteration 0 Tue Jul 8 00:30:19 2025
=-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul 8 00:30:21 2025
=-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul 8 00:30:23 2025
=-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul 8 00:30:26 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul 8 00:30:28 2025
=-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul 8 00:30:30 2025
=-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul 8 00:30:32 2025
=-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul 8 00:30:34 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul 8 00:30:37 2025
=-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul 8 00:30:39 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 1.6 4.8
Prediction (s.d.) 1.6 (0.36) 4.9 (1.07)
95% CI [0.95, 2.4] [2.84, 7.1]
Absolute effect (s.d.) -0.043 (0.36) -0.130 (1.07)
95% CI [-0.76, 0.65] [-2.28, 1.95]
Relative effect (s.d.) 3% (30%) 3% (30%)
95% CI [-32%, 69%] [-32%, 69%]
Posterior tail-area probability p: 0.44954
Posterior prob. of a causal effect: 55%
For more details, type: summary(impact, "report")## ── Estado: nivel local + estacionalidad + AR(10) ────────────────────────────
ss2d_qc <- list()
ss2d_qc <- AddLocalLevel(
ss2d_qc, c(as.numeric(df35_qc$y), rep(NA, 7))
)
ss2d_qc <- AddSeasonal(
ss2d_qc, c(as.numeric(df35_qc$y), rep(NA, 7)),
nseasons = 7, season.duration = 1
)
ss2d_qc <- AddAutoAr(
ss2d_qc, y = c(as.numeric(df35_qc$y), rep(NA, 7)),
lags = 10
)
## ── Covariables como data.frame ──────────────────────────────────────────────
covariates_qc <- df35_qc |> dplyr::select(
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin,
Y_lag10
)
covariates_future_qc <- df35_post_qc |> dplyr::select(
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin,
Y_lag10
)
covariates_full_qc <- dplyr::bind_rows(covariates_qc, covariates_future_qc)
## ── Regresión dinámica ───────────────────────────────────────────────────────
ss2d_qc <- AddDynamicRegression(ss2d_qc, covariates_full_qc)
## ── Ajuste BSTS y CausalImpact ───────────────────────────────────────────────
model2d_qc <- bsts(
c(as.numeric(df35_qc$y), rep(NA, 7)),
state.specification = ss2d_qc,
family = "student",
niter = 3e4,
seed = 2125
)
impact_qc <- CausalImpact(
bsts.model = model2d_qc,
model.args = list(prior.level.sd = .1, dynamic.regression = TRUE),
post.period.response = df35_post_qc$y
)
plot(impact_qc, "original")
summary(impact_qc)=-=-=-=-= Iteration 0 Tue Jul 8 00:30:45 2025
=-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul 8 00:30:48 2025
=-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul 8 00:30:51 2025
=-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul 8 00:30:54 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul 8 00:30:57 2025
=-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul 8 00:31:00 2025
=-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul 8 00:31:04 2025
=-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul 8 00:31:07 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul 8 00:31:10 2025
=-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul 8 00:31:13 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 1.4 10.0
Prediction (s.d.) 1.5 (0.46) 10.6 (3.23)
95% CI [0.62, 2.4] [4.35, 17.1]
Absolute effect (s.d.) -0.079 (0.46) -0.554 (3.23)
95% CI [-1.0, 0.81] [-7.1, 5.70]
Relative effect (s.d.) 6.2% (216%) 6.2% (216%)
95% CI [-42%, 126%] [-42%, 126%]
Posterior tail-area probability p: 0.43895
Posterior prob. of a causal effect: 56%
For more details, type: summary(impact, "report")## ── Vector de respuesta ───────────────────────────────────────────────────────
y35_qc <- c(as.numeric(df35_qc$y), rep(NA, 7))
## ── Covariables como data.frame ──────────────────────────────────────────────
covariates_qc <- df35_qc |> dplyr::select(
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin,
Y_lag10
)
covariates_future_qc <- df35_post_qc |> dplyr::select(
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin,
Y_lag10
)
covariates_full_qc <- dplyr::bind_rows(covariates_qc, covariates_future_qc)
## ── Especificación del estado: tendencia local lineal ────────────────────────
ss2e_qc <- list()
ss2e_qc <- AddLocalLinearTrend(ss2e_qc, y35_qc)
ss2e_qc <- AddSeasonal(ss2e_qc, y35_qc, nseasons = 7, season.duration = 1)
ss2e_qc <- AddAutoAr(ss2e_qc, y35_qc, lags = 1)
ss2e_qc <- AddDynamicRegression(ss2e_qc, covariates_full_qc)
## ── Ajuste BSTS y CausalImpact ───────────────────────────────────────────────
model35_qc <- bsts(
y35_qc,
state.specification = ss2e_qc,
family = "student",
niter = 3e4,
seed = 2125
)
impact2_qc <- CausalImpact(
bsts.model = model35_qc,
model.args = list(prior.level.sd = .1, dynamic.regression = TRUE),
post.period.response = df35_post_qc$y
)
summary(impact2_qc)
plot(impact2_qc)
=-=-=-=-= Iteration 0 Tue Jul 8 00:31:20 2025
=-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul 8 00:31:22 2025
=-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul 8 00:31:25 2025
=-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul 8 00:31:27 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul 8 00:31:30 2025
=-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul 8 00:31:32 2025
=-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul 8 00:31:35 2025
=-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul 8 00:31:37 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul 8 00:31:40 2025
=-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul 8 00:31:42 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 1.4 10.0
Prediction (s.d.) 1.5 (0.5) 10.8 (3.5)
95% CI [0.56, 2.6] [3.94, 17.9]
Absolute effect (s.d.) -0.11 (0.5) -0.76 (3.5)
95% CI [-1.1, 0.87] [-7.9, 6.10]
Relative effect (s.d.) -0.74% (1070%) -0.74% (1070%)
95% CI [-44%, 147%] [-44%, 147%]
Posterior tail-area probability p: 0.41321
Posterior prob. of a causal effect: 59%
For more details, type: summary(impact, "report")## ── Vector de respuesta ───────────────────────────────────────────────────────
y33_qc <- c(as.numeric(df33_qc$y), rep(NA, 7)) # 7 post-días + 2 de ajuste
## ── Covariables (data.frame) ─────────────────────────────────────────────────
covariates33_qc <- df33_qc |> dplyr::select(
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin,
Y_lag10
)
covariates_future33_qc <- df33_post_qc |> dplyr::select(
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin,
Y_lag10
)
covariates_full33_qc <- dplyr::bind_rows(covariates33_qc, covariates_future33_qc[1:7,])
## ── Especificación del estado: nivel local ───────────────────────────────────
ss2a_qc <- list()
ss2a_qc <- AddLocalLevel(ss2a_qc, y33_qc)
ss2a_qc <- AddSeasonal(ss2a_qc, y33_qc, nseasons = 7, season.duration = 1)
ss2a_qc <- AddAutoAr (ss2a_qc, y33_qc, lags = 10)
ss2a_qc <- AddDynamicRegression(ss2a_qc, covariates_full33_qc)
## ── Ajuste BSTS y CausalImpact ───────────────────────────────────────────────
model2a_qc <- bsts(
y33_qc,
state.specification = ss2a_qc,
family = "student",
niter = 3e4,
seed = 2125
)
impact33_qc <- CausalImpact(
bsts.model = model2a_qc,
model.args = list(prior.level.sd = .1, dynamic.regression = TRUE),
post.period.response = df33_post_qc$y[1:7]
)
plot(impact33_qc, "original")
summary(impact33_qc)=-=-=-=-= Iteration 0 Tue Jul 8 00:31:49 2025
=-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul 8 00:31:52 2025
=-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul 8 00:31:55 2025
=-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul 8 00:31:58 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul 8 00:32:01 2025
=-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul 8 00:32:04 2025
=-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul 8 00:32:07 2025
=-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul 8 00:32:10 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul 8 00:32:13 2025
=-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul 8 00:32:16 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 1.4 9.8
Prediction (s.d.) 1.5 (0.4) 10.2 (2.8)
95% CI [0.68, 2.3] [4.76, 15.8]
Absolute effect (s.d.) -0.066 (0.4) -0.460 (2.8)
95% CI [-0.86, 0.72] [-6.03, 5.03]
Relative effect (s.d.) -3.9% (1355%) -3.9% (1355%)
95% CI [-38%, 104%] [-38%, 104%]
Posterior tail-area probability p: 0.43406
Posterior prob. of a causal effect: 57%
For more details, type: summary(impact, "report")## ── Vector de respuesta ───────────────────────────────────────────────────────
y33_lin_qc <- c(as.numeric(df33_qc$y), rep(NA, 7))
## ── Covariables (data.frame) ─────────────────────────────────────────────────
covariates_full33_qc <- dplyr::bind_rows(covariates33_qc, covariates_future33_qc)
## ── Estado: tendencia local lineal ───────────────────────────────────────────
ss2b_qc <- list()
ss2b_qc <- AddLocalLinearTrend(ss2b_qc, y33_lin_qc)
ss2b_qc <- AddSeasonal (ss2b_qc, y33_lin_qc, nseasons = 7, season.duration = 1)
ss2b_qc <- AddAutoAr (ss2b_qc, y33_lin_qc, lags = 1)
ss2b_qc <- AddDynamicRegression(ss2b_qc, covariates_full33_qc)
## ── Ajuste BSTS y CausalImpact ───────────────────────────────────────────────
model33b_qc <- bsts(
y33_lin_qc,
state.specification = ss2b_qc,
family = "student",
niter = 3e4,
seed = 2125
)
impact2_33_qc <- CausalImpact(
bsts.model = model33b_qc,
model.args = list(prior.level.sd = .1, dynamic.regression = TRUE),
post.period.response = df33_post_qc$y[1:7]
)
summary(impact2_33_qc)
plot(impact2_33_qc)
=-=-=-=-= Iteration 0 Tue Jul 8 00:32:23 2025
=-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul 8 00:32:25 2025
=-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul 8 00:32:28 2025
=-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul 8 00:32:30 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul 8 00:32:32 2025
=-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul 8 00:32:35 2025
=-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul 8 00:32:37 2025
=-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul 8 00:32:39 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul 8 00:32:42 2025
=-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul 8 00:32:44 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 1.4 9.8
Prediction (s.d.) 1.5 (0.43) 10.7 (3.00)
95% CI [0.68, 2.4] [4.73, 16.7]
Absolute effect (s.d.) -0.12 (0.43) -0.86 (3.00)
95% CI [-0.98, 0.72] [-6.87, 5.06]
Relative effect (s.d.) 3.9% (164%) 3.9% (164%)
95% CI [-41%, 105%] [-41%, 105%]
Posterior tail-area probability p: 0.38691
Posterior prob. of a causal effect: 61%
For more details, type: summary(impact, "report")## ── Vector de respuesta ───────────────────────────────────────────────────────
y33b_qc <- c(as.numeric(df33_qc$y), rep(NA, 3)) # 7 post-días + 2 de ajuste
covariates_full33b_qc <- dplyr::bind_rows(covariates33_qc, covariates_future33_qc[1:3,])
## ── Especificación del estado: nivel local ───────────────────────────────────
ss2ab_qc <- list()
ss2ab_qc <- AddLocalLevel(ss2ab_qc, y33b_qc)
ss2ab_qc <- AddSeasonal(ss2ab_qc, y33b_qc, nseasons = 7, season.duration = 1)
ss2ab_qc <- AddAutoAr (ss2ab_qc, y33b_qc, lags = 10)
ss2ab_qc <- AddDynamicRegression(ss2ab_qc, covariates_full33b_qc)
## ── Ajuste BSTS y CausalImpact ───────────────────────────────────────────────
model2ab_qc <- bsts(
y33_qc,
state.specification = ss2ab_qc,
family = "student",
niter = 3e4,
seed = 2125
)Error in bsts(y33_qc, state.specification = ss2ab_qc, family = "student", : Caught exception with the following error message:
incompatible vector in SparseVector dot product:
dense vector: 0.0359117 3.55378e-06 -4.29923e-06 9.23099e-08 1.82078e-07 1.1973e-07 4.18066e-07 0 0 0 0 3.27532e-06 -6.80385e-06 4.98775e-06 -2.6448e-06 -1.53716e-06 9.6989e-07 -0.00153241 -0.000415958 -0.000115887 -0.000321038
sparse[0] = 1
sparse[1] = 1
sparse[7] = 1impact33b_qc <- CausalImpact(
bsts.model = model2ab_qc,
model.args = list(prior.level.sd = .1, dynamic.regression = TRUE),
post.period.response = df33_post_qc$y[1:3]
)Error: objeto 'model2ab_qc' no encontradoplot(impact33b_qc, "original")Error: objeto 'impact33b_qc' no encontradosummary(impact33b_qc)Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': objeto 'impact33b_qc' no encontrado## ── Vector de respuesta ───────────────────────────────────────────────────────
y33b_lin_qc <- c(as.numeric(df33_qc$y), rep(NA, 3))
## ── Covariables (data.frame) ─────────────────────────────────────────────────
covariates_full33b_qc <- dplyr::bind_rows(covariates33_qc, covariates_future33b_qc)Error: objeto 'covariates_future33b_qc' no encontrado## ── Estado: tendencia local lineal ───────────────────────────────────────────
ss2bb_qc <- list()
ss2bb_qc <- AddLocalLinearTrend(ss2bb_qc, y33b_lin_qc)
ss2bb_qc <- AddSeasonal (ss2bb_qc, y33b_lin_qc, nseasons = 7, season.duration = 1)
ss2bb_qc <- AddAutoAr (ss2bb_qc, y33b_lin_qc, lags = 1)
ss2bb_qc <- AddDynamicRegression(ss2bb_qc, covariates_full33b_qc)
## ── Ajuste BSTS y CausalImpact ───────────────────────────────────────────────
model33b_qc <- bsts(
y33b_lin_qc,
state.specification = ss2bb_qc,
family = "student",
niter = 3e4,
seed = 2125
)
impact2_33b_qc <- CausalImpact(
bsts.model = model33bb_qc,
model.args = list(prior.level.sd = .1, dynamic.regression = TRUE),
post.period.response = df33_post_qc$y[1:3]
)Error: objeto 'model33bb_qc' no encontradosummary(impact2_33b_qc)Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': objeto 'impact2_33b_qc' no encontradoplot(impact2_33b_qc)Error: objeto 'impact2_33b_qc' no encontrado=-=-=-=-= Iteration 0 Tue Jul 8 00:32:52 2025
=-=-=-=-=
=-=-=-=-= Iteration 3000 Tue Jul 8 00:32:54 2025
=-=-=-=-=
=-=-=-=-= Iteration 6000 Tue Jul 8 00:32:56 2025
=-=-=-=-=
=-=-=-=-= Iteration 9000 Tue Jul 8 00:32:58 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Jul 8 00:33:00 2025
=-=-=-=-=
=-=-=-=-= Iteration 15000 Tue Jul 8 00:33:03 2025
=-=-=-=-=
=-=-=-=-= Iteration 18000 Tue Jul 8 00:33:05 2025
=-=-=-=-=
=-=-=-=-= Iteration 21000 Tue Jul 8 00:33:07 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Jul 8 00:33:09 2025
=-=-=-=-=
=-=-=-=-= Iteration 27000 Tue Jul 8 00:33:11 2025
=-=-=-=-=cat(paste0("R library: ", Sys.getenv("R_LIBS_USER")))
cat(paste0("Date: ",withr::with_locale(new = c('LC_TIME' = 'C'), code =Sys.time())))
cat(paste0("Editor context: ", getwd()))
cat("quarto version: "); system("quarto --version")
quarto::quarto_version()
save.image("_data/step3.RData")R library: H:/Mi unidad/PERSONAL ANDRES/UCH_salud_publica/pasantia/f1/f1/renv/library/windows/R-4.4/x86_64-w64-mingw32Date: 2025-07-08 00:33:14.740894Editor context: H:/Mi unidad/PERSONAL ANDRES/UCH_salud_publica/pasantia/f1/f1quarto version: [1] 0
[1] '1.6.39'sesion_info <- devtools::session_info()Warning in system2("quarto", "-V", stdout = TRUE, env = paste0("TMPDIR=", : el
comando ejecutado '"quarto"
TMPDIR=C:/Users/andre/AppData/Local/Temp/RtmpYrN5xg/file569457bd5dcb -V' tiene
el estatus 1dplyr::select(
tibble::as_tibble(sesion_info$packages),
c(package, loadedversion, source)
) |>
knitr::kable(caption = "R packages", format = "html",
col.names = c("Row number", "Package", "Version"),
row.names = FALSE,
align = c("c", "l", "r")) |>
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12)|>
kableExtra::scroll_box(width = "100%", height = "375px") | Row number | Package | Version |
|---|---|---|
| abind | 1.4-8 | RSPM |
| assertthat | 0.2.1 | RSPM |
| backports | 1.5.0 | RSPM |
| bayesplot | 1.12.0 | RSPM |
| bdsmatrix | 1.3-7 | RSPM |
| bit | 4.6.0 | CRAN (R 4.4.3) |
| bit64 | 4.6.0-1 | CRAN (R 4.4.3) |
| Boom | 0.9.15 | RSPM |
| BoomSpikeSlab | 1.2.6 | RSPM |
| boot | 1.3-30 | CRAN (R 4.4.1) |
| bpmn | 0.1.0 | Github (bergant/bpmn@628d3efaa27544c221b2fd7f1895301a63b70c49) |
| bridgesampling | 1.1-2 | RSPM |
| brms | 2.22.0 | RSPM |
| Brobdingnag | 1.2-9 | RSPM |
| broom | 1.0.8 | RSPM |
| bsts | 0.9.10 | RSPM |
| cachem | 1.1.0 | CRAN (R 4.4.3) |
| car | 3.1-3 | RSPM |
| carData | 3.0-5 | RSPM |
| CausalImpact | 1.3.0 | RSPM |
| CBPS | 0.23 | RSPM |
| checkmate | 2.3.2 | RSPM |
| cli | 3.6.5 | RSPM |
| coda | 0.19-4.1 | RSPM |
| codetools | 0.2-20 | CRAN (R 4.4.1) |
| collapse | 2.1.2 | RSPM |
| colorspace | 2.1-1 | RSPM |
| conquer | 1.3.3 | RSPM |
| cubature | 2.1.4 | RSPM |
| curl | 6.2.3 | CRAN (R 4.4.1) |
| CVXR | 1.0-15 | RSPM |
| dagitty | 0.3-4 | RSPM |
| data.table | 1.17.4 | RSPM |
| devtools | 2.4.5 | RSPM |
| DHARMa | 0.4.7 | RSPM |
| digest | 0.6.37 | RSPM |
| distributional | 0.5.0 | RSPM |
| doParallel | 1.0.17 | RSPM |
| doRNG | 1.8.6.2 | RSPM |
| doSNOW | 1.0.20 | RSPM |
| dplyr | 1.1.4 | RSPM |
| dreamerr | 1.5.0 | RSPM |
| ECOSolveR | 0.5.5 | RSPM |
| ellipsis | 0.3.2 | RSPM |
| emmeans | 1.11.1 | RSPM |
| estimability | 1.5.1 | RSPM |
| evaluate | 1.0.3 | RSPM |
| farver | 2.1.2 | RSPM |
| fastDummies | 1.7.5 | RSPM |
| fastmap | 1.2.0 | CRAN (R 4.4.3) |
| fect | 1.0.0 | RSPM |
| fixest | 0.12.1 | RSPM |
| forcats | 1.0.0 | RSPM |
| foreach | 1.5.2 | RSPM |
| forecast | 8.24.0 | RSPM |
| Formula | 1.2-5 | RSPM |
| fracdiff | 1.5-3 | RSPM |
| fs | 1.6.6 | RSPM |
| future | 1.49.0 | RSPM |
| future.apply | 1.11.3 | RSPM |
| geeM | 0.10.1 | RSPM |
| geepack | 1.3.12 | RSPM |
| generics | 0.1.4 | RSPM |
| geosphere | 1.5-20 | RSPM |
| GGally | 2.2.1 | RSPM |
| ggdag | 0.2.13 | RSPM |
| ggplot2 | 3.5.2 | RSPM |
| ggstats | 0.9.0 | RSPM |
| glmmTMB | 1.1.11 | RSPM |
| glmnet | 4.1-9 | RSPM |
| glmx | 0.2-1 | RSPM |
| globals | 0.18.0 | RSPM |
| glue | 1.8.0 | RSPM |
| gmp | 0.7-5 | RSPM |
| gnm | 1.1-5 | RSPM |
| gridExtra | 2.3 | RSPM |
| gtable | 0.3.6 | RSPM |
| gtools | 3.9.5 | RSPM |
| hms | 1.1.3 | CRAN (R 4.4.3) |
| htmltools | 0.5.8.1 | RSPM |
| htmlwidgets | 1.6.4 | RSPM |
| httpuv | 1.6.16 | RSPM |
| httr2 | 1.1.2 | RSPM |
| igraph | 2.1.4 | RSPM |
| iterators | 1.0.14 | RSPM |
| jsonlite | 2.0.0 | CRAN (R 4.4.3) |
| kableExtra | 1.4.0 | RSPM |
| kernlab | 0.9-33 | RSPM |
| knitr | 1.50 | RSPM |
| labeling | 0.4.3 | RSPM |
| later | 1.4.2 | RSPM |
| lattice | 0.22-6 | CRAN (R 4.4.1) |
| lfe | 3.1.1 | RSPM |
| lifecycle | 1.0.4 | RSPM |
| listenv | 0.9.1 | RSPM |
| lme4 | 1.1-37 | RSPM |
| lmtest | 0.9-40 | RSPM |
| loo | 2.8.0 | RSPM |
| lubridate | 1.9.4 | RSPM |
| magrittr | 2.0.3 | RSPM |
| MASS | 7.3-60.2 | CRAN (R 4.4.1) |
| MatchIt | 4.7.2 | RSPM |
| Matrix | 1.7-0 | CRAN (R 4.4.1) |
| MatrixModels | 0.5-4 | RSPM |
| matrixStats | 1.5.0 | RSPM |
| maxLik | 1.5-2.1 | RSPM |
| memoise | 2.0.1 | CRAN (R 4.4.3) |
| mgcv | 1.9-1 | CRAN (R 4.4.1) |
| mime | 0.13 | CRAN (R 4.4.3) |
| miniUI | 0.1.2 | RSPM |
| minqa | 1.2.8 | RSPM |
| miscTools | 0.6-28 | RSPM |
| mvtnorm | 1.3-3 | RSPM |
| nixtlar | 0.6.2 | RSPM |
| nlme | 3.1-164 | CRAN (R 4.4.1) |
| nloptr | 2.2.1 | RSPM |
| nnet | 7.3-19 | CRAN (R 4.4.1) |
| np | 0.60-18 | RSPM |
| numDeriv | 2016.8-1.1 | RSPM |
| openxlsx | 4.2.8 | RSPM |
| optimx | 2025-4.9 | RSPM |
| PanelMatch | 3.1.1 | RSPM |
| parallelly | 1.44.0 | RSPM |
| pillar | 1.10.2 | RSPM |
| pkgbuild | 1.4.8 | RSPM |
| pkgconfig | 2.0.3 | RSPM |
| pkgload | 1.4.0 | RSPM |
| plm | 2.6-6 | RSPM |
| plyr | 1.8.9 | RSPM |
| posterior | 1.6.1 | RSPM |
| pracma | 2.4.4 | CRAN (R 4.4.1) |
| processx | 3.8.6 | RSPM |
| profvis | 0.4.0 | RSPM |
| promises | 1.3.2 | RSPM |
| ps | 1.9.1 | RSPM |
| purrr | 1.0.4 | RSPM |
| Qtools | 1.5.9 | RSPM |
| quadprog | 1.5-8 | RSPM |
| quantdr | 1.2.2 | RSPM |
| quantmod | 0.4.28 | RSPM |
| quantreg | 6.1 | RSPM |
| quarto | 1.4.4 | RSPM |
| qvcalc | 1.0.4 | RSPM |
| R6 | 2.6.1 | RSPM |
| rappdirs | 0.3.3 | CRAN (R 4.4.3) |
| rbibutils | 2.3 | RSPM |
| RColorBrewer | 1.1-3 | RSPM |
| Rcpp | 1.0.14 | RSPM |
| RcppParallel | 5.1.10 | RSPM |
| Rdpack | 2.6.4 | RSPM |
| readr | 2.1.5 | CRAN (R 4.4.3) |
| reformulas | 0.4.1 | RSPM |
| relimp | 1.0-5 | RSPM |
| remotes | 2.5.0 | RSPM |
| renv | 1.1.2 | CRAN (R 4.4.1) |
| reshape2 | 1.4.4 | RSPM |
| rgenoud | 5.9-0.11 | RSPM |
| rio | 1.2.3 | RSPM |
| rlang | 1.1.6 | RSPM |
| rmarkdown | 2.29 | RSPM |
| Rmpfr | 1.1-0 | RSPM |
| rngtools | 1.5.2 | RSPM |
| rstantools | 2.4.0 | RSPM |
| rstudioapi | 0.17.1 | RSPM |
| sandwich | 3.1-1 | RSPM |
| scales | 1.4.0 | RSPM |
| scpi | 3.0.0 | RSPM |
| sessioninfo | 1.2.3 | RSPM |
| shape | 1.4.6.1 | RSPM |
| shiny | 1.10.0 | RSPM |
| snow | 0.4-4 | RSPM |
| sp | 2.2-0 | RSPM |
| SparseM | 1.84-2 | RSPM |
| stringi | 1.8.7 | RSPM |
| stringmagic | 1.2.0 | RSPM |
| stringr | 1.5.1 | RSPM |
| survival | 3.6-4 | CRAN (R 4.4.1) |
| svglite | 2.2.1 | RSPM |
| Synth | 1.1-8 | RSPM |
| systemfonts | 1.2.3 | RSPM |
| tensorA | 0.36.2.1 | RSPM |
| textshaping | 1.0.1 | RSPM |
| tibble | 3.2.1 | RSPM |
| tidygraph | 1.3.1 | RSPM |
| tidyr | 1.3.1 | RSPM |
| tidyselect | 1.2.1 | RSPM |
| tidyverse | 2.0.0 | RSPM |
| timechange | 0.3.0 | RSPM |
| timeDate | 4041.110 | RSPM |
| TMB | 1.9.17 | RSPM |
| tseries | 0.10-58 | RSPM |
| TTR | 0.24.4 | RSPM |
| tzdb | 0.5.0 | CRAN (R 4.4.3) |
| urca | 1.3-4 | RSPM |
| urlchecker | 1.0.1 | RSPM |
| usethis | 3.1.0 | RSPM |
| V8 | 6.0.3 | RSPM |
| vctrs | 0.6.5 | RSPM |
| viridisLite | 0.4.2 | RSPM |
| weathercan | 0.7.3.9000 | https://ropensci.r-universe.dev (R 4.4.3) |
| withr | 3.0.2 | RSPM |
| xfun | 0.52 | RSPM |
| xml2 | 1.3.8 | CRAN (R 4.4.3) |
| xtable | 1.8-4 | RSPM |
| xts | 0.14.1 | RSPM |
| yaml | 2.3.10 | RSPM |
| zip | 2.3.3 | RSPM |
| zoo | 1.8-14 | RSPM |
reticulate::py_list_packages()%>%
knitr::kable(caption = "Python packages", format = "html",
col.names = c("Package", "Version", "Requirement"),
row.names = FALSE,
align = c("c", "l", "r", "r"))%>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12)|>
kableExtra::scroll_box(width = "100%", height = "375px") | Package | Version | Requirement |
|---|---|---|
| absl-py | 2.1.0 | absl-py==2.1.0 |
| asttokens | 2.4.1 | asttokens==2.4.1 |
| astunparse | 1.6.3 | astunparse==1.6.3 |
| audioconverter | 2.0.3 | audioconverter==2.0.3 |
| autograd | 1.6.2 | autograd==1.6.2 |
| autograd-gamma | 0.5.0 | autograd-gamma==0.5.0 |
| beautifulsoup4 | 4.12.3 | beautifulsoup4==4.12.3 |
| Brotli | 1.1.0 | Brotli==1.1.0 |
| certifi | 2023.11.17 | certifi==2023.11.17 |
| cffi | 1.16.0 | cffi==1.16.0 |
| charset-normalizer | 3.3.2 | charset-normalizer==3.3.2 |
| clarabel | 0.9.0 | clarabel==0.9.0 |
| click | 8.1.7 | click==8.1.7 |
| cloudpickle | 3.0.0 | cloudpickle==3.0.0 |
| colorama | 0.4.6 | colorama==0.4.6 |
| comm | 0.2.1 | comm==0.2.1 |
| contourpy | 1.2.0 | contourpy==1.2.0 |
| cvxopt | 1.3.2 | cvxopt==1.3.2 |
| cvxpy | 1.5.2 | cvxpy==1.5.2 |
| cycler | 0.12.1 | cycler==0.12.1 |
| debugpy | 1.8.0 | debugpy==1.8.0 |
| decorator | 4.4.2 | decorator==4.4.2 |
| delete-chrome-history-py | 0.1.8 | delete-chrome-history-py==0.1.8 |
| easyocr | 1.7.1 | easyocr==1.7.1 |
| ecos | 2.0.13 | ecos==2.0.13 |
| editdistance | 0.8.1 | editdistance==0.8.1 |
| efficientnet | 1.0.0 | efficientnet==1.0.0 |
| essential-generators | 1.0 | essential-generators==1.0 |
| et-xmlfile | 1.1.0 | et-xmlfile==1.1.0 |
| executing | 2.0.1 | executing==2.0.1 |
| fancyimpute | 0.7.0 | fancyimpute==0.7.0 |
| ffmpeg | 1.4 | ffmpeg==1.4 |
| ffmpeg-python | 0.2.0 | ffmpeg-python==0.2.0 |
| filedir | 0.0.3 | filedir==0.0.3 |
| filelock | 3.13.1 | filelock==3.13.1 |
| flatbuffers | 24.3.25 | flatbuffers==24.3.25 |
| fonttools | 4.47.2 | fonttools==4.47.2 |
| formulaic | 1.0.1 | formulaic==1.0.1 |
| fsspec | 2023.12.2 | fsspec==2023.12.2 |
| future | 0.18.3 | future==0.18.3 |
| gast | 0.6.0 | gast==0.6.0 |
| git-filter-repo | 2.45.0 | git-filter-repo==2.45.0 |
| google-pasta | 0.2.0 | google-pasta==0.2.0 |
| graphviz | 0.20.3 | graphviz==0.20.3 |
| grpcio | 1.65.4 | grpcio==1.65.4 |
| gTTS | 2.5.1 | gTTS==2.5.1 |
| h5py | 3.11.0 | h5py==3.11.0 |
| idna | 3.6 | idna==3.6 |
| imageio | 2.34.2 | imageio==2.34.2 |
| imageio-ffmpeg | 0.5.1 | imageio-ffmpeg==0.5.1 |
| imgaug | 0.4.0 | imgaug==0.4.0 |
| iniconfig | 2.0.0 | iniconfig==2.0.0 |
| interface-meta | 1.3.0 | interface-meta==1.3.0 |
| ipykernel | 6.29.5 | ipykernel==6.29.5 |
| ipython | 8.20.0 | ipython==8.20.0 |
| jedi | 0.19.1 | jedi==0.19.1 |
| Jinja2 | 3.1.3 | Jinja2==3.1.3 |
| joblib | 1.4.0 | joblib==1.4.0 |
| jupyter_client | 8.6.0 | jupyter_client==8.6.0 |
| jupyter_core | 5.7.1 | jupyter_core==5.7.1 |
| keras | 3.4.1 | keras==3.4.1 |
| Keras-Applications | 1.0.8 | Keras-Applications==1.0.8 |
| keras-ocr | 0.9.3 | keras-ocr==0.9.3 |
| kiwisolver | 1.4.5 | kiwisolver==1.4.5 |
| knnimpute | 0.1.0 | knnimpute==0.1.0 |
| lazy_loader | 0.4 | lazy_loader==0.4 |
| libclang | 18.1.1 | libclang==18.1.1 |
| lifelines | 0.28.0 | lifelines==0.28.0 |
| llvmlite | 0.41.1 | llvmlite==0.41.1 |
| Markdown | 3.6 | Markdown==3.6 |
| markdown-it-py | 3.0.0 | markdown-it-py==3.0.0 |
| MarkupSafe | 2.1.4 | MarkupSafe==2.1.4 |
| matplotlib | 3.8.2 | matplotlib==3.8.2 |
| matplotlib-inline | 0.1.6 | matplotlib-inline==0.1.6 |
| mdurl | 0.1.2 | mdurl==0.1.2 |
| mido | 1.3.3 | mido==1.3.3 |
| ml-dtypes | 0.4.0 | ml-dtypes==0.4.0 |
| more-itertools | 10.2.0 | more-itertools==10.2.0 |
| moviepy | 1.0.3 | moviepy==1.0.3 |
| mpmath | 1.3.0 | mpmath==1.3.0 |
| multipledispatch | 1.0.0 | multipledispatch==1.0.0 |
| mutagen | 1.47.0 | mutagen==1.47.0 |
| namex | 0.0.8 | namex==0.0.8 |
| natsort | 8.4.0 | natsort==8.4.0 |
| nest-asyncio | 1.5.9 | nest-asyncio==1.5.9 |
| networkx | 3.2.1 | networkx==3.2.1 |
| ninja | 1.11.1.1 | ninja==1.11.1.1 |
| nose | 1.3.7 | nose==1.3.7 |
| numba | 0.58.1 | numba==0.58.1 |
| numexpr | 2.10.0 | numexpr==2.10.0 |
| numpy | 1.26.3 | numpy==1.26.3 |
| openai-whisper | 20231117 | openai-whisper==20231117 |
| opencv-python | 4.10.0.84 | opencv-python==4.10.0.84 |
| opencv-python-headless | 4.10.0.84 | opencv-python-headless==4.10.0.84 |
| openpyxl | 3.1.4 | openpyxl==3.1.4 |
| opt-einsum | 3.3.0 | opt-einsum==3.3.0 |
| optree | 0.12.1 | optree==0.12.1 |
| osqp | 0.6.5 | osqp==0.6.5 |
| packaging | 23.2 | packaging==23.2 |
| pandas | 2.2.0 | pandas==2.2.0 |
| pandas-flavor | 0.6.0 | pandas-flavor==0.6.0 |
| parso | 0.8.3 | parso==0.8.3 |
| patsy | 0.5.6 | patsy==0.5.6 |
| pillow | 10.2.0 | pillow==10.2.0 |
| platformdirs | 4.1.0 | platformdirs==4.1.0 |
| pluggy | 1.5.0 | pluggy==1.5.0 |
| polars | 1.9.0 | polars==1.9.0 |
| proglog | 0.1.10 | proglog==0.1.10 |
| prompt-toolkit | 3.0.43 | prompt-toolkit==3.0.43 |
| protobuf | 4.25.4 | protobuf==4.25.4 |
| psutil | 5.9.8 | psutil==5.9.8 |
| pure-eval | 0.2.2 | pure-eval==0.2.2 |
| pyarrow | 15.0.0 | pyarrow==15.0.0 |
| pyclipper | 1.3.0.post5 | pyclipper==1.3.0.post5 |
| pycparser | 2.22 | pycparser==2.22 |
| pycryptodomex | 3.20.0 | pycryptodomex==3.20.0 |
| pydotplus | 2.0.2 | pydotplus==2.0.2 |
| pydub | 0.24.1 | pydub==0.24.1 |
| Pygments | 2.17.2 | Pygments==2.17.2 |
| pyjanitor | 0.26.0 | pyjanitor==0.26.0 |
| PyMuPDF | 1.24.9 | PyMuPDF==1.24.9 |
| PyMuPDFb | 1.24.9 | PyMuPDFb==1.24.9 |
| pyparsing | 3.1.1 | pyparsing==3.1.1 |
| PyPDF2 | 3.0.1 | PyPDF2==3.0.1 |
| pyreadr | 0.5.0 | pyreadr==0.5.0 |
| pytesseract | 0.3.10 | pytesseract==0.3.10 |
| pytest | 8.3.1 | pytest==8.3.1 |
| python-bidi | 0.6.0 | python-bidi==0.6.0 |
| python-dateutil | 2.8.2 | python-dateutil==2.8.2 |
| pytube | 15.0.0 | pytube==15.0.0 |
| pytube3 | 9.6.4 | pytube3==9.6.4 |
| pytz | 2023.3.post1 | pytz==2023.3.post1 |
| pywin32 | 306 | pywin32==306 |
| PyYAML | 6.0.1 | PyYAML==6.0.1 |
| pyzmq | 25.1.2 | pyzmq==25.1.2 |
| qdldl | 0.1.7.post1 | qdldl==0.1.7.post1 |
| regex | 2023.12.25 | regex==2023.12.25 |
| requests | 2.32.3 | requests==2.32.3 |
| rich | 13.7.1 | rich==13.7.1 |
| rpy2 | 3.5.16 | rpy2==3.5.16 |
| scikit-image | 0.24.0 | scikit-image==0.24.0 |
| scikit-learn | 1.3.2 | scikit-learn==1.3.2 |
| scikit-survival | 0.22.2 | scikit-survival==0.22.2 |
| scipy | 1.11.4 | scipy==1.11.4 |
| scs | 3.2.6 | scs==3.2.6 |
| seaborn | 0.13.2 | seaborn==0.13.2 |
| semantic-version | 2.10.0 | semantic-version==2.10.0 |
| setuptools-rust | 1.8.1 | setuptools-rust==1.8.1 |
| shapely | 2.0.5 | shapely==2.0.5 |
| six | 1.16.0 | six==1.16.0 |
| soupsieve | 2.5 | soupsieve==2.5 |
| SpeechRecognition | 3.10.1 | SpeechRecognition==3.10.1 |
| spyder-kernels | 2.5.2 | spyder-kernels==2.5.2 |
| stack-data | 0.6.3 | stack-data==0.6.3 |
| statsmodels | 0.14.1 | statsmodels==0.14.1 |
| sympy | 1.12 | sympy==1.12 |
| target | 0.0.11 | target==0.0.11 |
| tensorboard | 2.17.0 | tensorboard==2.17.0 |
| tensorboard-data-server | 0.7.2 | tensorboard-data-server==0.7.2 |
| tensorflow | 2.17.0 | tensorflow==2.17.0 |
| tensorflow-intel | 2.17.0 | tensorflow-intel==2.17.0 |
| tensorflow-io-gcs-filesystem | 0.31.0 | tensorflow-io-gcs-filesystem==0.31.0 |
| termcolor | 2.4.0 | termcolor==2.4.0 |
| threadpoolctl | 3.4.0 | threadpoolctl==3.4.0 |
| tifffile | 2024.7.24 | tifffile==2024.7.24 |
| tiktoken | 0.5.2 | tiktoken==0.5.2 |
| torch | 2.4.0 | torch==2.4.0 |
| torchaudio | 2.4.0 | torchaudio==2.4.0 |
| torchvision | 0.19.0 | torchvision==0.19.0 |
| tornado | 6.4 | tornado==6.4 |
| tqdm | 4.66.1 | tqdm==4.66.1 |
| traitlets | 5.14.1 | traitlets==5.14.1 |
| translator | 0.0.9 | translator==0.0.9 |
| typing_extensions | 4.9.0 | typing_extensions==4.9.0 |
| tzdata | 2023.4 | tzdata==2023.4 |
| tzlocal | 5.2 | tzlocal==5.2 |
| urllib3 | 2.1.0 | urllib3==2.1.0 |
| validators | 0.33.0 | validators==0.33.0 |
| watchdog | 3.0.0 | watchdog==3.0.0 |
| wcwidth | 0.2.13 | wcwidth==0.2.13 |
| websockets | 12.0 | websockets==12.0 |
| Werkzeug | 3.0.3 | Werkzeug==3.0.3 |
| whisper | 1.1.10 | whisper==1.1.10 |
| wrapt | 1.16.0 | wrapt==1.16.0 |
| xarray | 2024.1.1 | xarray==2024.1.1 |
| youtube-dl | 2021.12.17 | youtube-dl==2021.12.17 |
| yt-dlp | 2024.7.9 | yt-dlp==2024.7.9 |